18LX

Wind trajectory

Show code
library(GeoPressureR)
library(tidyverse)
library(leaflet)
library(leaflet.extras)
library(raster)
library(dplyr)
library(ggplot2)
library(plotly)
knitr::opts_chunk$set(echo = FALSE)
load(paste0("../data/1_pressure//", params$gdl_id, "_pressure_prob.Rdata"))
load(paste0("../data/2_light/", params$gdl_id, "_light_prob.Rdata"))
load(paste0("../data/3_static/", params$gdl_id, "_static_prob.Rdata"))
# load(paste0("../data/4_basic_graph/", params$gdl_id, "_basic_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_wind_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_grl.Rdata"))
col <- rep(RColorBrewer::brewer.pal(8, "Dark2"), times = ceiling(max(tag$sta$sta_id) / 8))

Altitude

Altitudes are computed based on pressure measurement of the geolocation, corrected based on the assumed location of the shortest path. This correction accounts therefore for the natural variation of pressure as estimated by ERA-5. The vertical lines indicate the sunrise (dashed) and sunset (solid).

Show code
p <- ggplot() +
  # geom_line(data = tag$pressure, aes(x = date, y = obs), colour = "grey") +
  geom_line(data = do.call("rbind", shortest_path_timeserie), aes(x = date, y = altitude)) +
  geom_line(data = do.call("rbind", shortest_path_timeserie) %>% filter(sta_id > 0), aes(x = date, y = altitude, col = factor(sta_id))) +
  geom_vline(data = twl, aes(xintercept = twilight, linetype = ifelse(rise, "dashed", "solid"), color="grey"), lwd=0.1) +
  theme_bw() +
  scale_colour_manual(values = col) +
  scale_y_continuous(name = "Altitude (m)")

ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)

Wintering location

Show code
file <- paste0("figure_print/wintering_location/wintering_location_",params$gdl_id,".png")
if(file.exists(file)){
  knitr::include_graphics(file)
}

Latitude time

Show code
 tmp <- lapply(pressure_prob, function(x) {
    mt <- metadata(x)
    df <- data.frame(
      start = mt$temporal_extent[1],
      end = mt$temporal_extent[2],
      sta_id = mt$sta_id
    )
  })
  tmp2 <- do.call("rbind", tmp)

sim_lat <- as.data.frame(t(path_sim$lat)) %>%
  mutate(sta_id = path_sim$sta_id) %>%
  pivot_longer(-c(sta_id)) %>%
  left_join(tmp2,by="sta_id")

sim_lat_p <- sim_lat %>%
  filter(sta_id==max(sta_id)) %>%
  mutate(start=end) %>%
  rbind(sim_lat)

sp_lat <- as.data.frame(shortest_path) %>% left_join(tmp2,by="sta_id")

sp_lat_p <- sp_lat %>%
  filter(sta_id==max(sta_id)) %>%
  mutate(start=end) %>%
  rbind(sp_lat)

p <- ggplot() +
  geom_step(data=sim_lat_p, aes(x=start, y=value, group=name), alpha=.07) +
  geom_point(data=sp_lat_p, aes(x=start, y=lat)) +
  xlab('Date') +
  ylab('Latitude') +
  theme_light()

ggplotly(p, dynamicTicks = T)

Shortest path and simulated path

The large circles indicates the shortest path (overall most likely trajectory) estimated by the graph approach. The size is proportional to the duration of stay. The small dots and grey lines represents 10 possible trajeectories of the bird according to the model.

Click on the full-screen mode button on the top-left of the map to see more details on the map.

Show code
sta_duration <- unlist(lapply(static_prob_marginal, function(x) {
  as.numeric(difftime(metadata(x)$temporal_extent[2], metadata(x)$temporal_extent[1], units = "days"))
}))
pal <- colorFactor(col, as.factor(seq_len(length(col))))
m <- leaflet(width = "100%") %>%
  addProviderTiles(providers$Stamen.TerrainBackground) %>%
  addFullscreenControl() %>%
  addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
  addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = pal(factor(shortest_path$sta_id, levels = tag$sta$sta_id)), weight = sta_duration^(0.3) * 10)

for (i in seq_len(nrow(path_sim$lon))) {
  m <- m %>%
    addPolylines(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = 0.5, weight = 1, color = "#808080") %>%
    addCircles(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = .7, weight = 1, color = pal(factor(shortest_path$sta_id, levels = tag$sta$sta_id)))
}
m

Marginal probability map

The marginal probability map estimate the overall probability of position at each stationary period regardless of the trajectory taken by the bird. It is the most useful quantification of the uncertainty of the position of the bird.

Show code
li_s <- list()
l <- leaflet(width = "100%") %>%
  addProviderTiles(providers$Stamen.TerrainBackground) %>%
  addFullscreenControl()
for (i_r in seq_len(length(static_prob_marginal))) {
  i_s <- metadata(static_prob_marginal[[i_r]])$sta_id
  info <- metadata(static_prob_marginal[[i_r]])$temporal_extent
  info_str <- paste0(i_s, " | ", info[1], "->", info[2])
  li_s <- append(li_s, info_str)
  l <- l %>%
    addRasterImage(static_prob_marginal[[i_r]], colors = "OrRd", opacity = 0.8, group = info_str) %>%
    addCircles(lng = shortest_path$lon[i_s], lat = shortest_path$lat[i_s], opacity = 1, color = "#000", weight = 10, group = info_str)
}
l %>%
  addLayersControl(
    overlayGroups = li_s,
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  hideGroup(tail(li_s, length(li_s) - 1))

Wind assistance

Show code
  fun_marker_color <- function(norm){
    if (norm < 20){
      "darkpurple"
    } else if (norm < 35){
      "darkblue"
    } else if (norm < 50){
      "lightblue"
    } else if (norm < 60){
      "lightgreen"
    } else if (norm < 80){
      "yellow"
    } else if (norm < 100){
      "lightred"
    } else {
      "darkred"
    }
  }
  fun_NSEW <- function(angle){
    angle <- angle  %% (pi* 2)
    angle <- angle*180/pi
    if (angle < 45/2){
      "E"
    } else if (angle < 45*3/2){
      "NE"
    } else if (angle < 45*5/2){
      "N"
    } else if (angle < 45*7/2){
      "NW"
    } else if (angle < 45*9/2){
      "W"
    } else if (angle < 45*11/2){
      "SW"
    } else if (angle < 45*13/2){
      "S"
    }else if (angle < 45*15/2){
      "SE"
    } else {
      "E"
    }
  }

  sta_duration <- unlist(lapply(static_prob_marginal,function(x){as.numeric(difftime(metadata(x)$temporal_extent[2],metadata(x)$temporal_extent[1],units="days"))}))

  m <-leaflet(width = "100%") %>%
    addProviderTiles(providers$Stamen.TerrainBackground) %>%  addFullscreenControl() %>%
    addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
    addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#000", weight = sta_duration^(0.3)*10)

  for (i_s in seq_len(grl$sz[3]-1)){
    if (grl$flight_duration[i_s]>5){
      edge <- which(grl$s == shortest_path$id[i_s] & grl$t == shortest_path$id[i_s+1])

      label = paste0( i_s,': ', grl$flight[[i_s]]$start, " - ", grl$flight[[i_s]]$end, "<br>",
                      "F. dur.: ", round(grl$flight_duration[i_s]), ' h <br>',
                      "GS: ", round(abs(grl$gs[edge])), ' km/h, ',fun_NSEW(Arg(grl$gs[edge])),'<br>',
                      "WS: ", round(abs(grl$ws[edge])), ' km/h, ',fun_NSEW(Arg(grl$ws[edge])),'<br>',
                      "AS: ", round(abs(grl$as[edge])), ' km/h, ',fun_NSEW(Arg(grl$as[edge])),'<br>')

      iconArrow <- makeAwesomeIcon(icon = "arrow-up",
                                   library = "fa",
                                   iconColor = "#FFF",
                                   iconRotate = (90 - Arg(grl$ws[edge])/pi*180) %% 360,
                                   squareMarker = TRUE,
                                   markerColor = fun_marker_color(abs(grl$ws[edge])))

      m <- m %>% addAwesomeMarkers(lng = (shortest_path$lon[i_s] + shortest_path$lon[i_s+1])/2,
                                   lat = (shortest_path$lat[i_s] + shortest_path$lat[i_s+1])/2,
                                   icon = iconArrow, popup = label)
    }
  }
  m

Histogram of Speed

Show code
edge <- t(graph_path2edge(path_sim$id, grl))
nj <- ncol(edge)
nsta <- ncol(path_sim$lon)

speed_df <- data.frame(
  as = abs(grl$as[edge]),
  gs = abs(grl$gs[edge]),
  ws = abs(grl$ws[edge]),
  sta_id_s = rep(head(grl$sta_id,-1), nj),
  sta_id_t = rep(tail(grl$sta_id,-1), nj),
  flight_duration = rep(head(grl$flight_duration,-1), nj),
  dist = geosphere::distGeo(
    cbind(as.vector(t(path_sim$lon[,1:nsta-1])), as.vector(t(path_sim$lat[,1:nsta-1]))),
    cbind(as.vector(t(path_sim$lon[,2:nsta])),   as.vector(t(path_sim$lat[,2:nsta])))
  ) / 1000
) %>% mutate(
  name = paste(sta_id_s,sta_id_t, sep="-")
)

plot1 <- ggplot(speed_df, aes(reorder(name, sta_id_s), gs)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot2 <- ggplot(speed_df, aes(reorder(name, sta_id_s), ws)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot3 <- ggplot(speed_df, aes(reorder(name, sta_id_s), as)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot4 <- ggplot(speed_df, aes(reorder(name, sta_id_s), flight_duration)) + geom_point() + theme_bw() +scale_x_discrete(name = "")

subplot(ggplotly(plot1), ggplotly(plot2), ggplotly(plot3), ggplotly(plot4), nrows=4, titleY=T)

Table of transition

Show code
alt_df = do.call("rbind", shortest_path_timeserie) %>%
    arrange(date) %>%
    mutate(
      sta_id_s = cummax(sta_id),
      sta_id_t = sta_id_s+1
    ) %>%
    filter(sta_id == 0 & sta_id_s > 0 ) %>%
    group_by(sta_id_s, sta_id_t) %>%
    summarise(
      alt_min = min(altitude),
      alt_max = max(altitude),
      alt_mean = mean(altitude),
      alt_med = median(altitude),
    )

  trans_df <- speed_df  %>%
    group_by(sta_id_s,sta_id_t,flight_duration) %>%
    summarise(
      as_m = mean(as),
      as_s = sd(as),
      gs_m = mean(gs),
      gs_s = sd(gs),
      ws_m = mean(ws),
      ws_s = sd(ws),
      dist_m = mean(dist),
      dist_s = sd(dist)
    ) %>%
    left_join(alt_df)

trans_df %>% kable()
sta_id_s sta_id_t flight_duration as_m as_s gs_m gs_s ws_m ws_s dist_m dist_s alt_min alt_max alt_mean alt_med
1 2 3.4166667 34.57011 17.494468 47.68620 21.711983 18.003929 3.8248656 162.92785 74.18261 157.25331 1060.0105 654.1848 604.5139
2 3 7.3333333 37.18249 15.983395 44.39281 16.646386 11.408376 5.5214162 325.54730 122.07350 240.07618 985.7225 595.2773 547.5521
3 4 7.9166667 40.58991 15.977018 56.06959 18.303820 19.826997 5.5771796 443.88426 144.90524 -230.29992 1404.8480 298.0325 214.8520
4 5 4.9166667 33.81206 17.370380 34.80925 16.913821 4.123158 2.5650512 171.14546 83.15962 22.06228 1241.4214 487.8901 506.4481
5 6 9.8333333 41.63452 11.463145 50.63115 16.145250 11.780202 2.6895206 497.87295 158.76162 9.47344 1670.3913 679.8558 587.9694
6 7 9.4166667 36.63912 13.033330 36.98587 12.525579 12.530809 5.0572737 348.28359 117.94921 127.30825 3458.0496 1838.6165 1684.6414
7 9 28.4166667 55.26767 17.586456 64.81036 14.110518 11.531872 3.7032344 1841.69430 400.97388 NA NA NA NA
9 10 8.5000000 62.58095 19.842185 71.09937 18.510293 16.076736 2.4089642 604.34461 157.33749 935.15268 2205.3068 1804.8507 1854.0280
10 11 8.9166667 41.04600 17.261752 48.03535 18.698191 20.129450 4.1477964 428.31520 166.72553 658.94857 3100.9877 1691.4406 1707.8348
11 12 5.2500000 41.59498 21.239332 40.42210 19.901930 3.833050 1.4844871 212.21603 104.48513 295.45738 1527.9584 1134.8679 1336.0644
12 13 7.1666667 44.52462 17.710289 38.32131 17.694810 8.304144 1.8732288 274.63602 126.81281 359.33122 1215.7742 706.7753 650.7422
13 14 0.8333333 35.96519 16.702519 36.62206 20.024884 9.981134 3.0560772 30.51839 16.68740 136.05294 219.7510 177.9019 177.9019
14 16 12.6666667 12.42522 6.712675 12.08427 7.529162 3.545367 0.9657721 153.06746 95.36939 NA NA NA NA
16 17 3.7500000 32.95460 11.883280 32.75673 11.478415 2.103072 1.1820068 122.83775 43.04406 220.08353 411.3244 317.3691 321.5518
17 18 1.5000000 33.64262 10.438046 39.00593 17.020175 17.288656 1.6683386 58.50889 25.53026 141.33298 171.3591 154.2831 150.1572
18 19 10.6666667 37.47178 11.796753 43.69214 10.826953 18.739869 2.9062068 466.04952 115.48749 243.59053 1672.0438 896.6158 895.9195
19 20 4.9166667 33.69398 14.300953 42.22981 17.339492 19.528778 4.9322794 207.62988 85.25250 312.52187 1782.4563 906.9879 825.7417
20 21 4.2500000 36.99632 14.036391 28.68552 18.504806 25.468923 4.2639659 121.91347 78.64543 287.69319 1124.9316 655.0057 669.5111
21 22 10.6666667 36.38271 17.266830 30.44620 12.582570 16.071009 7.0997736 324.75946 134.21408 432.44483 1905.1577 794.3203 628.2020
22 23 10.5000000 39.79834 16.936986 38.94978 16.247836 12.807723 5.7257082 408.97265 170.60228 789.05124 3154.8776 2129.4470 1842.2562
23 24 10.2500000 39.65864 13.748037 78.32582 15.589634 43.528736 5.2998547 802.83962 159.79375 1035.30185 3786.9226 3129.5720 3261.6373
24 26 23.5000000 51.53299 8.173692 113.85681 3.927738 66.987229 6.5230024 2675.63511 92.30185 NA NA NA NA
26 27 7.9166667 26.75544 7.335790 31.91815 8.167588 8.769154 1.2363085 252.68532 64.66007 340.90962 1229.2793 706.1691 737.0251
27 28 5.7500000 39.82505 16.393176 57.14521 21.375869 20.570207 6.4380061 328.58494 122.91125 474.38752 2804.7508 1480.3557 1389.4821
28 29 6.5000000 61.13500 15.352408 62.17359 17.253713 25.230900 5.6896699 404.12837 112.14913 207.11052 2572.7330 964.6875 450.2296